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In this paper we propose a method for the construction of locally conservative flux 

fields from Generalized Multiscale Finite Element Method (GMsFEM) pressure 

solutions. The flux values are obtained from an element-based postprocessing 

procedure in which an independent set of 4 x 4 linear systems need to be solved. 

2 To test the performance of the method we consider two heterogeneous permeabil- 

£j ity coefficients and couple the resulting fluxes to a two-phase flow model. The 

increase in accuracy associated with the computation of the GMsFEM pressure 

^ solutions is inherited by the postprocessed flux fields and saturation solutions, 

CN and is closely correlated to the size of the reduced-order systems. In particular, 

^ the addition of more basis functions to the enriched coarse space yields solutions 

\0 that more accurately capture the behavior of the fine scale model. A number of 

'^ numerical examples are offered to validate the performance of the method. 
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Many physical processes in science and engineering are described by partial 
differential equations whose coefficients vary over many length scales. Typical 
examples may include subsurface flows where the permeability of the porous 
medium is represented by a high-contrast, heterogeneous coefficient. In recent 
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decades, multiscale methods have been introduced as an effective tool for treat- 
ing these types of problems [1, 2, 12, 16, 20, 21]. An important component of 
this class of methods is the independent construction of a set of multiscale basis 
functions that span a solution space that is tied to a coarse grid, i.e., one whose 
discretization parameter is much larger than the characteristic scale of the hetero- 
geneous coefficient. In particular, once a precomputed set of basis functions is 
available, a specified global coupling mechanism may be used in order to obtain 
the associated coarse scale solution. As the fine scale information is imbedded 
into the basis functions, a coarse grid solution inherits the fine scale effects of the 
underlying system. In other words, the multiscale basis functions offer a direct 
method of projecting a coarse solution to the fine grid. The present paper con- 
siders a class of multiscale methods that will be used to effectively solve elliptic 
pressure equations that appear in a two-phase flow model. 

While standard multiscale methods have proven effective for a variety of ap- 
plications (see, e.g., [11, 16, 21]), we employ a more recent framework in which 
the coarse space may be systematically enriched so that the approximate solu- 
tion sought in it converges to the fine grid solution. The enrichment procedure 
hinges on the construction of localized spectral problems, where dominant eigen- 
functions are used in the construction of the enriched space [10, 15]. This type 
of spectral enrichment allows for the number basis functions (and the size of the 
coarse space) to be flexibly chosen such that a desired level of numerical accuracy 
may be achieved. This framework, which is coined as the Generalized Multiscale 
Finite Element Method (GMsFEM), incorporates the enriched solution space into 
a continuous Galerkin (CG) global formulation in order to obtain approximate 
pressure solutions. 

An advantage of employing a continuous Galerkin multiscale formulation is 
the relative ease of implementation and resemblence to standard finite element 
variational formulation. However, a well known limitation of CG is that the re- 
sulting solution does not satisfy local conservation. In particular, in the cases 
when it is necessary to couple the resulting fluxes to a transport equation, lo- 
cal conservation is required. While finite volume-type methods, mixed methods, 
and discontinuous Galerkin methods typically guarantee conservation [1, 7, 11], 
the respective formulations yield systems that are more delicate (and sometimes 
larger) than the CG counterpart. Furthermore, to the best of our knowledge the 
blend of enrichment techniques with multiscale methods are still at its infancy as 
there has not been any attempt to carry out the formulation using other than contin- 
uous Galerkin formulation (see, e.g., [9] for a recent development using discontin- 
uous Galerkin method). As a result, we consider the alternative of postprocessing 



a global CG solution in order to obtain the desired conservation. Numerous meth- 
ods have been proposed in order to postprocess finite element solutions to obtain 
conservative fluxes (see, e.g., [6, 19, 22, 24]), however, in this paper we general- 
ize the procedure offered in [4] in which the authors perform a global solve and 
subsequent element-based computations to achieve conservation. 

In this paper we propose a technique that provides flux conservation in the 
context of GMsFEM. For rectangular finite elements, our method hinges on a 
postprocessing technique in which independent 4x4 systems of equations are 
solved on each coarse element to obtain the conservative fluxes. We note that 
similar derivation can be accomplished for triangular finite elements that yields 
an independent 3x3 system of equations. While the postprocessing procedure 
yields conservative fluxes on the coarse scale (which might suffice for some target 
applications), we also employ an independent downscaling procedure to construct 
a conservative flux field on the underlying fine grid. We note that coarse scale con- 
servative discontinuous Galerkin GMsFEM formulations have been used in (see, 
e.g., [9]), yet emphasize that the method proposed in this paper requires no mod- 
ification to the original CG formulation and allows for the fluxes to be computed 
on the fine grid. Furthermore, to our knowledge, conservative GMsFEM-type 
methods have not yet been incorporated for solving multi-phase flow models in 
the existing literature. To test the performance of the proposed method we solve 
a standard two-phase flow model using distinct cases of high-contrast permeabil- 
ity coefficients. In all cases, an increase in the dimension of the coarse solution 
space yields solutions that are shown to more accurately capture the behavior of 
the fine scale. In particular, the error decline of the elliptic solution (which has 
been rigorously analyzed in [10]), is directly inherited by the resulting flux values 
and two-phase saturation solutions. 

The rest of the paper is organized as follows. In Sect. 2 we introduce a standard 
two-phase model along with a description of the operator splitting technique that 
is used for solving the model. In Sect. 3 we describe the Generalized Multiscale 
Finite Element Method (GMsFEM), and follow the construction by introducing 
the procedure for the computation of postprocessed, conservative flux quantities 
in Sect. 4. A variety of numerical tests are offered in Sect. 5 in order to validate the 
performance the proposed method. To finish the paper we offer some concluding 
remarks in Sect. 6. 



2. Model problem 

2.7. Two-phase model 

We consider a heterogeneous oil reservoir which is confined in a domain Q. 
The reservoir is equipped with an injection well, from which water is discharged 
to displace the trapped oil towards the production wells, situated on the perimeter 
of the domain. The dynamics of the movement of the fluids in the reservoir are 
categorized as an immiscible two-phase system with water and oil (denoted by w 
and o, respectively) that is incompressible. Capillary pressure is not included in 
the model. Further simplifying assumptions that we use are a gravity-free environ- 
ment and that the two fluids fill the pore space. Then, the Darcy's law combined 
with a statement of conservation of mass allow us to write the governing equations 
of the flow as 

V • V = ^, where v = -A{S)k{xy^p, (2.1) 

and 

— + V.(/(5)v) = ^,, (2.2) 

where v is the Darcy velocity, S is the water saturation, and k is the permeability 
coefficient. The total mobility A{S) and the flux function f{S) are respectively 
given by: 

i/c\ ^rw\S) kro\S) krw{S)/fiw .^ ^. 

^(^)= + . JK^) = 7777T . (2.3) 

where krj, j = w, o, is the relative permeability of the phase j. 

2.2. Solution algorithm 

Notice that the elliptic part (2.1) and the transport part (2.2) of the system are 
coupled through the total mobility. In order to solve this problem numerically, 
we use an operator splitting technique [3], where saturation at the previous time 
step is used when solving the elliptic part of the system to obtain the velocity v. 
This velocity is then used in an explicit time stepping scheme for the transport 
equation. This velocity is held constant for a predetermined number of time steps, 
which yields a new saturation. This new saturation is then used to solve the elliptic 
problem again and the process is continued until the final time is reached. A 
schematic of the operator splitting is shown in Fig. 1. 

To discretize the transport equation, we first integrate (2.2) with respect to 
time and over some Q c Q. Here we apply the left end point quadrature rule to 
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Figure 1 : Operator splitting for the two-phase problem. 

its second term in time and use integration by parts to arrive at the approximation 

meas(C,)(5,,^-5,,^_i) + A^ I v • w/(5,,^_i)d/ = I q^d.x, 

Jdc, Jc, 

where we have neglected the error terms and 

1 r 

^ z,n ~ 77TT I ^ v^' ^n) dX. 

meas(C,) Jq 

The transport part of the system is solved with an explicit time stepping as seen 
above and an upwinding scheme is used on the term J v • nf(S ) d/. A review of 
up winding on a rectangular mesh can be found for example in [25]. Typically in 
this situation, it is imperative that numerical approximation of v satisfies certain 
local conservation property. In particular, given V/^ ^ v, it is desirable to have 

I Vh'ndl= j qdx. (2.4) 

A natural way of obtaining this local conservation property is to seek the solution 
of (/?, v) from the first order system (2.1). At the approximation stage, one ends 
up with the mixed finite element formulation [23]. A more common approach is 
to transform (2.1) into a second order equation that governs p. The approximate 
solution for p is sought after which the approximate solution of v is calculated 
using the relation v = -A(S)k(x)Vp, and hence a postprocessing procedure is 
used. Unfortunately, a standard technique such as the Galerkin finite element 
method does not allow for a straightforward postprocessing to obtain a Vh that is 
locally conservative. 

Furthermore, aside from the issues pertaining to local conservation, it is almost 
always impossible to conduct numerical simulations at the finer scales if one is to 
include the ever increasing details of geological information. A viable alternative 
is the use of multiscale methods, which is the subject of the next section. 



3. Generalized multiscale finite element method 

In this section we describe a systematic coarse grid solution technique that 
may be used as a reduced-order alternative to a standard fine grid approach (such 
as a fully resolved finite element discretization). The solution procedure is built 
within the framework of the Generalized Multiscale Finite Element Method (GMs- 
FEM), in which a pressure solution is sought within a space of precomputed mul- 
tiscale basis functions. This type of generalized method allows for the systematic 
enrichment of coarse solution spaces based on the underlying structure of the 
problem (e.g., the structure of the permeability coefficient k(x)). 

To fix our attention, we consider 

f -V • (Ak(x)Vp) = q in Q c M^ 

p = Pd on Fd (3.1) 

-Ak(x)Vp ' n = g^ on Fn 

with Dirichlet and Neumann boundary conditions given by Pb^En^ respectively, 
and a forcing function q. Within the context of the operator splitting procedure 
we assume that A is already available. The variational formulation associated with 
(3.1) is to seek p with (p - p^y) e H^^ = {v e H^(Q) : vIfd = 0} that satisfies 

a(p, w) = ((?, w) + (gN, w)rN "^weH^, (3.2) 

where 

a(p,w) = I Ak(x)Vp 'Vwdx, (q,w) = I qwdx, and (gN.>^)rN = I gN>^d/. 
Jq Jq Jfn 

Typical Galerkin finite element methods seek the approximate solution of p 
in some finite dimensional subspace of H^ that satisfies (3.2). This finite di- 
mensional subspace is associated with a discretization of Q into Th consisting 
of closed quadrilateral (or triangular) elements r such that Q = U^e^^r, where 
h = max^e^J/^^} and hj is the diameter of r. For example, by defining the con- 
forming linear finite element space "Vh over Th as 

"^h = \^h ^ C(Q) : Whir is linear for all r eTh and Whlvj^ = Ol, 

the approximate solution ph is found to satisfy (ph - PB,h) ^ "^h and 

ci{ph. ^h) = (<?, ^h) + <gN, w/,)rN V w/, G ^h' (3.3) 



To proceed, we designate Z as the set of all vertices in the discretization of Q 
where Z = Zin U Z^ U Z^ with Z^ is the set of vertices on Fd, Zin is the set of 
interior vertices, and Z^ is the set of vertices on Fn- Designating {^^l^eZinUZn ^s the 
linear/bilinear basis functions ofVh^ (3.3) yields a linear system 

Ap = f, (3.4) 

where A is a square matrix with entries a(0^, 0^), f is a vector with entries (/, 0^) 
for ^ G Zin and (/, 0^) + (g^,^^) for ^ e Z^. The vector p contains the nodal 
values of /7/^, i.e., /7^ = Ph(x^), which coincide with the coordinates in the linear 
combination of {0^}. We note that the size of this system is dim(Zin U Zn). 

3.1. MsFEM for pressure equations 

The multiscale finite element method (MsFEM) was introduced in [17] and 
further analyzed in [18]. Similar to the approximation method described earlier, 
MsFEM is a continuous Galerkin finite element method that is based on solving 
(3.2) in a finite dimensional subspace of H^. Its distinct feature is in the careful 
choice of a multiscale finite dimensional subspace that allows calculation of the 
approximate solution on Th without directly resolving the fine scale heterogeneity 
globally. This means h is much larger than the characteristic fine scale of k(x). 
Since much of the fine scale heterogeneity has its source from k(x), this informa- 
tion should be ingrained in this subspace. This is done through the construction 
of the multiscale basis functions that characterize the subspace. 

As in the standard Galerkin finite element method, the multiscale basis func- 
tions are associated with the nodes/vertices in the discretization of Q. In the inter- 
est of off'ering a straightforward presentation, in what follows we assume that Th 
is a collection of rectangles r such as depicted in Fig. 2. Extension of the method 
to triangles follows the same line of arguments. 

For a vertex z g Zin U Zn, the corresponding multiscale basis functions ;t/^ are 
defined in such a way ihai Xz,t - Xz\ is governed by 

-V . {k{x)Vx,,r) = 0, in T 

XzA^) = 0z,t(^), on 5t, 

where 02,t(x) is the standard bilinear function in r and z is a vertex of r. The 
multiscale finite dimensional space is defined as 

^n.s,/z = spanj;^, : z g Zin U Zn) c //^. (3.6) 



The MsFEM solution is /^ms,/? with (pms,h - PD,h) ^ '^ms,h satisfying 

Ci(Pms,h. W/z) = (<?, W/z) + (gN, W/,)rN V W/, G ^ms,h' (3.7) 

Similar to the standard Galerkin finite element method, the linear system resulting 
from (3.7) is 

AmsP = fms, (3.8) 

where A^s is a square matrix with entries a(x^,Xz)^ fms is a vector with entries 
(f.Xd fo^ f ^ ^in and (/,/\f^) + (g^^x^) for f G Zn, and p is as in (3.4). 

To this end, we emphasize on the significant advantage of MsFEM. If one is 
to solve the pressure equation (3.1) to the extent that the fine scale heterogeneity 
of k(x) is directly resolved in the finite element formulation (3.3), then ideally the 
level of mesh resolution on which the flow and transport are simulated has to be 
in a comparable scale with the level of subsurface resolution exhibited by k(x). In 
turn, the resulting linear system such as expressed in (3.4) can have a very large 
dimension whose inversion poses a very challenging task. Employing MsFEM on 
the other hand, avoids this drawback as there is no need to pose (3.7) on a mesh 
having comparable size to the characteristic scale of k(x). This is because the 
multiscale finite element space in (3.6) already contains this information. 

Still, potentially there can be a significant error associated with MsFEM ap- 
proximation that stems from the imposition of linear boundary condition onxz on 
dr, which causes a mismatch as compared to the true solution. This mismatch 
is even more pronounced when k(x) exhibits channelized features and/or scat- 
tered inclusions in which the the values of k(x) are orders of magnitude higher 
(or lower) compared to the neighboring regions. Within the context of modeling 
and simulation of multiphase flow and transport, there have been several research 
directions aiming toward alleviating this situation. For example, oversampling 
techniques [5], adoption of global information [11], and the use of a local-global 
iterative approach [8] have been used for error reduction. More recent investiga- 
tions involve the alternative of systematically enriching the original coarse space 
n^nis,/z- This is the subject of next subsection. 

3.2. GMsFEM for pressure equations 

The Generalized Multiscale Finite Element Method (GMsFEM) is based on a 
systematic enrichment of H^ms,/? [10, 14, 15]. This enrichment is made available by 
taking advantage of the knowledge of the spectral properties of the original differ- 
ential operator governing the multiscale basis functions ;i;^; namely, the left hand 
side of (3.5). These types of enriched spaces yield pressure solutions whose errors 

8 



decrease with respect to the locahzed eigenvalue behavior. In related work, it is 
shown that the errors typically depend on the first eigenvalue that is not included 
in the space construction. We refer the interested reader to [10] for rigorous error 
estimates. 
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Figure 2: Discretization ofQ. into Th = Ur. Here oj^ = u1 ^r^ is the suppO^^) 



To equip the description below, for a vertex z ^ Zin U Zn, we let oj^ be the 
support of the multiscale basis function ;i;^, namely, oj^ = ^%i^j^ where Tj are 
finite elements having z as one of its vertices. Enrichment of "Vms,/? employs the 
pointwise energy of the original multiscale basis functions 



k = k J] h'wxf. 



ZGZinUZn 



(3.9) 



In particular, using k as data, we solve an eigenvalue problem 

I -V • (fc(x)V(Aj = ju,kil/„ in CO, 
I -VjAz ■ n = 0, on doj,. 



(3.10) 



for each a>,. We denote the eigenvalues and eigenvectors of (3.10) by {/^^/j and 
{jA^/}, respectively. By direct observation of (3.10) we see that the first eigenpair 
is /^2 1 = and jAz,i = 1- We order the resulting eigenvalues as 



jtl^j < ytl^,2 < . . . < yUz,„ < . . . 



(3.11) 



The size of the eigenvalues is closely related to the structure of k and, in particu- 
lar, m inclusions and channels in k yields m asymptotically vanishing eigenvalues. 
We use the eigenvectors corresponding to small, asymptotically vanishing eigen- 
values for the construction of the enriched space. In particular, we define the basis 
functions 

®z/ = Xz^z,i for z G Zin U Zn and \<i< L„ (3.12) 

where L^ denotes the number of eigenvectors that will be chosen for each vertex 
z. We note that this setting yields supp(0^/) = co^. With the updated multiscale 
basis functions available, we define the enriched multiscale finite element space 
as 

%„,s,/z = span{o,,, : z g Zin U Zn and 1 < ^ < L,) c Z/^. (3.13) 

The GMsFEM solution is /7ems,/z with (;7ems,/z - /^d,/z) ^ %ms,/z satisfying 

aiPomsM "^h) = (<?, Wh) + <gN, H^/z)rN V W/, G %ms,/z. (3.14) 

Upon comparison, it is obvious that "Vms,/? ^ 'T^ems,/z and thus dim(n^ems,/z) > 
dim('Vms,h)' Consequently, the resulting linear system AgmsP = fems has a larger di- 
mension than its counterpart in (3.8), where here we view p as a vector containing 
the coordinates (rather than the pressure nodal values) in the linear combination of 
O^ / that span the approximate pressure. However, this system is still significantly 
smaller compared to the fine scale analogue in (3.4) that is solved on mesh that has 
size of the order of the characteristic scale of ^(x). Thus, the enriched coarse sys- 
tem offers a suitable reduced-order alternative for obtaining approximate pressure 
solutions while maintaining an acceptable level of accuracy. 

We reiterate that the construction of the enriched basis functions in Eq. (3.12) 
is performed on a respective oj^ (refer back to Fig. 2). However, the postprocessing 
technique offered in the next section is localized to r. As such, it is important to 
note that the enriched basis functions need to be restricted into respective element 
contributions to implement the proposed procedure. So while we refer to the 
enriched basis functions synonymously (whether they are posed on a oj^ or r), we 
make this distinction for additional clarity. 

4. Locally conservative flux by postprocessing of the GMsFEM solution 

A pivotal contribution of this section is the introduction of a procedure in 
which enriched pressure solutions may be postprocessed to ensure that the con- 
servation property in (2.4) is met. In doing so, we may use the conservative flux 

10 



quantities that are required to solve the two-phase model described in Sect. 2. 
Additionally, we remark that the conservative flux quantities and associated satu- 
ration profiles are shown to inherit the increased level of numerical accuracy that 
is associated with the enriched coarse space construction. 

The postprocessing technique that we propose in order to construct?/^ • n sat- 
isfying (2.4) is achieved by relegating the evaluation to independent element-by- 
element calculations. At this stage, we need to make precise about what Q should 
look like. A suitable choice within the configuration of the nodally based Galerkin 
finite element method is to set Q as the control volume associated with vertex z 
(see left plot of Fig. 3). In this respect, we employ a basic fact about GMsFEM 
that we obtain from (3.14): 



a{ptmsM ^z,i) = (<?. ^z/). for Z G Zin, ^ = 1, • • • , L^ 



(4.1) 



Since supp(0^/) = C0z = ^ =i^ 



J' 



(^,0,,^)= ^(x)0,,^dx = J]f,,^j, (4.2) 

r ^ 

\JTf^ndio^ ■_j 



with 



Fz,ej= I q{x)^^,tdx, and G^^tj = j ^N<l>^/d/. 



(4.3) 



We conclude for each £ = 1, • • • , L^ that 

4 4 

X(ez/j-/^z/j) = 0, if zeZi„, J](Q,,ij-F,jj-G,,ej) = 0, if z € Z„. (4.4) 



7=1 7=1 
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Figure 3: Left: Q is the control volume associated with vertex z, where dC^ = E^^^ E^yU E^^^ E^^. 
Right: A finite element r is divided into four quadrilaterals t^,tx,ty, and t^. 



4.1. Auxiliary boundary value problem in r 

Since we desire independent element-by-element calculations, we proceed 
with a formulation of an auxiliary boundary value problem in r (see the right 
plots of Fig. 3) having vertices v(t) = {w, x,_y,z}. The boundary value problem 
governs pv with 

-V • {Xk{x)Vpr) = q in T 

-Ak(x)V^j ' n =^j on dr. 

Here we designate dr = U^evCr)^/. where E^. = dr n dt^ (i.e. half of each element 
edge containing the vertex ^). Furthermore, we set g^ as piecewise function on dr 
such that 

§; d/ = F^,i - 2^,1, for ^ G v(t). 



(4.5) 



I 



(4.6) 



where 



(4.7) 



2^,1 = I AkVp^rns,h ' VO^j dx and F^,i = I (?0^,i dx. 
The existence and uniqueness of the above problem is stated below 

Lemma 1. The compatibility condition holds for (4.5) and thus p is unique up to 
a constant. 

Proof. For the compatibility condition [13], it suffices to show that 

I grd/= I ^dx. 

Jdr Jt 
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Because ^ 0^,i| = ^ x^,t = 1, we get 

Xi ^^'1 " N Zi ^^'1 ^^ " I ^^-^ 

J] 2^,1 = r^/:V/7e„,s,/. • V( J^ 0^,i)dx = 0. 

Taking into consideration the choice of g^ above, these then give 

r ?r d/ = J] r ?r d/ = J] (F^,i - e^,i) = \q dx, (4.8) 

which is the desired result. D 

We note that (4.8) imphes 

- r AkVpr 'ndl = - f AkVp . fid/, (4.9) 

which shows that the solution of (4.5) recovers the flux ofp averaged over dr, i.e., 
a local conservation property. We use (4.5) as a governing principle to derive the 
postprocessing technique for calculating a locally conservative flux from /^ems,/?- 

4.2. Elemental calculation 

This elemental calculation is based on discretization of r into quadrilaterals t^, 
i.e., T = U^ev(T)^^, ^ach of which yields f ^ = Q Pi r, see the right plots of Fig. 3. 
We set the local solution space as "Vir) = span{0^ iJ^evCr)- The numerical solution 
associated with (4.5) is to find pj^h ^ "^i^) satisfying 



- I XkVpjh' ndl = I qdx, "i^evir). 

Jdtr Jtr 



(4.10) 



The following four equations result from (4.10): 



?wx - <fwz = 2w,l - -Fw.l + ? dx, 



I 



-^wx + ?Iy = Gx,l - i^x,l + I ^ dx, 

Jty 

<l\y^(fwz = Qz,i-F,,i + I qdx. 
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(4.11) 



where 



JEI^ J Ely 

ql, = - I AkVpr,h'ndl, qly = - I AkVpr,h'ndl, 

and E], = dt^ Pi dtj^, for ^, 77 = w, x, y, z, and ^ ^ //. We note that similar equations 
will be created if r is triangle (see [4]). 

Because pv,/? = Z^evCr) ^^^^,1 with unknown a^, the above equation yields a 
linear system of the form Aa = f where 

A^r^ = - I /IfcVO;^,! • fid/ and 7^= I (?dx- I g;d/. (4.12) 

Here the system is of dimension 4. We note that in the case of r adjacent to Fn, a 
similar linear system can be derived that takes into account the effect of g^. 

Since this linear system is obtained from numerical approximation of (4.5), 
which is a Neumann problem, the matrix A is singular. On the other hand, because 
the system has a small dimension, we may specify the value at one of the vertices 
to remove Jhe null space of A, for instance by adding a constant to one of the 
entries of A, and then invert the modified matrix. The fact that a is not unique is 
irrelevant because the desired final result from the system is the flux as governed 
by q^. which is unique. 

4.3. Upscaled local conservation 

The next lemma verifies that the aggregation of elemental calculations satisfies 
(2.4) for z e Zin. We note that the same is true for z g Zn whose proof we omit for 
simplicity. 

Lemma 2. Fix a z ^ Z/„ with co^ = u1 ^r^ and control volume Q (see Fig. 3). Set 
I Jh'ndl = q^^l + q^^l, I Vh' ndl = ql^ + q^ 

I Jh'ndl = ql\-\-ql\, I 7h ' ndl = ql' -\- ql 



T-zy 



fi- 



Then 



I Vk' ndl = I qdx. 
Jdc, Jc, 
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Proof. The basis of the proof is using the last equation in (4.11) applied to r^, for 
j = 1, • • • ,4. We perform a direct calculation and rearrange the terms involved to 
get 

I V/^ • fi d/ = I V/^ • fi d/ + I V/^ • fi d/ + I V/^ • fi d/ + I V/^ • n d/ 



4 4 ^ 



where we have appropriately translated the local indexing in (4.11) for Q^i and 
F^j into Q^^ij and /^^,ij, respectively. Notice that Zj=i(Gz,ij ~ ^z,ij) = by (4.4). 
Furthermore, t^j = Q n r^ and thus u1 ^f^ ^ = Q. This completes the proof. D 

As we see from the above description, what the postprocessing has been able 
to accomplish is a reconstruction of upscaled (or averaged) conservative flux in 
Q whose diameter is on the order of h. However, MsFEM (and for that matter 
GMsFEM) always utilizes h that is far greater than the characteristic scale of k{x). 
In this context, the upscaled flux gathered in Lemma 2 is comparable to entries 
of p coming from solving (3.7) (respectively from solving (3.14) for GMsFEM). 
Since "Vms,/? and "Vems,/? are built from the multiscale basis functions containing fine 
scale information of k{x), having this p allows for calculation of the approximate 
pressure down to the level of characteristic scale of k{x). However, the current 
stage of development does not yet allow the upscaled flux to have this capability. 

At many practical levels, obtaining upscaled locally conservative fluxes al- 
ready allows for simulation of multiphase flow and the transport problem in Sect. 2. 
However, the transport equation (2.2) must then be discretized at the same mesh 
level as where the approximate pressure is solved, i.e., with h that is far greater 
than the characteristic scale of k{x). Indeed, this practice is sufficient and suitable 
for some target applications. Nonetheless, for a large class of problems, such 
as modeling flows in channelized subsurface with potential localized features, 
achieving an acceptable accuracy does require the simulation to be performed 
with a discretization parameter that is comparable with the scale of k{x). For this 
to occur, we need the flux to be conservative at that finer scale. In the next sub- 
section, we offer a downscaling procedure that allows this capability by taking 
advantage of the upscaled conservative flux above. 
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4.4. A downscaling procedure 

The main driving fact for the downscahng procedure is the reahzation that 
after the postprocessing in Sect. 4, we have 



dc, Jc, 



which can be thought of a statement of as compatibihty condition in Q. This 
means we can proceed with formulating a boundary problem 

-V • (Ak(x)Vp^ ) = q in C. 

1^1 (4.14) 

-Ak(x)Vp^^ ' n = Vh' n on 5Q. 

Here v'h = -^k^PT,h that is evaluated pointwise on segments of dC^ that are in r. 
In fact, this is the same calculation that we did in order to derive (4.11). In this 
way the compatibility condition is readily satisfied, and thus the solution of (4.14) 
exists. Similar to the postprocessing in Sect. 4, nonuniqueness of the solution is 
of no concern since our interest is to gather -Ak^p^ from (4.14). 

In our implementation, we numerically solve (4.14) for every Q in Q with the 
discretization of Q using the same mesh configuration as that of r when numer- 
ically solving (3.5). Obviously, the associated mesh parameter should be smaller 
than h and comparable to the characteristic scale of k{x). We use the standard 
Galerkin finite element method (3.3) for both of these problems. In addition, the 
numerical solution of (4.14) is further postprocessed (see [4] for the description) 
to get a locally conservative flux on the finer scale. Once this is done for all Q, 
we obtain a the downscaled locally conservative flux in Q, which is in turn used 
in the simulation of transport equation (2.2). 

We should like to emphasize that the main advantage of the proposed series 
of upscaled and downscaled postprocessings is due to their independence of each 
other. The procedure is immediately parallelizable, fits well in the framework of 
CPU-GPU clusters, while the cost for communication is minimal. In this respect, 
they are indeed in the same spirit as the multiscale basis functions calculations. 

5. Numerical examples 

A variety of numerical examples that validate the performance of the proposed 
method are presented in this section. In particular, we perform a convergence 
study to address the convergence of the flux to a fine-scale reference solution as 
the number of enriched basis functions is increased. In general, the improvement 
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Figure 4: Permeabilities used in the examples. The plot on the left is shown in physical scale, while 
the right is shown in log scale with the left being referred to as deterministic (with /:max/^min ~ 
2 X 10"^) and the right as random (with ^max/^min ~ 1.8 x 10^). 



of the solution is shown to be dependent on the type of permeabihty coefficient 
and the level of enrichment. Applications to single and multi-phase flow are also 
presented in this section. 

We consider two distinct permeabilities within this section. Their spatial pro- 
files are shown in Fig. 4 with the left field shown in physical scale, and the right 
in log scale. All of the applications presented in this section use the domain 
Q = [0, 1] X [0, 1]. The left plot is a deterministic, high-contrast coefficient with 
abrupt transitions between regions of low and high permeability. This perme- 
ability is posed on a fine mesh of 100 x 100 elements. The right plot in Fig. 4 
is a single realization of a random, channelized permeability that is posed on a 
120 X 120 mesh. Both examples of permeability exhibit high-contrast features, 
which can make solving Eq. (2.1) a demanding task. The ratio between the max- 
imum (/:max) and minimum value (kmin) of k(x) can be thought of as representing 
the condition number of the resulting linear system in the finite element method. 
Compared with other methods such as the finite volume element method or mixed 
finite element method, continuous Galerkin finite element method has an advan- 
tage when combined with the postprocessing because the resulting linear systems 
can be easier to solve [4]. 

5.1. Single-phase flow 

In this subsection, we compare the velocities obtained from postprocessing 
Eq. (2.1) to be coupled to the transport problem in Eq. (2.2). To solve the pressure 
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Table 1 : Comparison of relative velocity errors for different levels of enrichment. The L^ norm 
of the downscaled velocity was computed for each of the cases below and compared with the 
fully-resolved reference values. 





Deterministic 


Random 


4 = 1 


1.458 


0.401 


4 = 2 


0.480 


0.320 


4 = 4 


0.203 


0.283 


4 = 6 


0.199 


0.274 



equation in (2.1) we use Dirichlet conditions /?l = 1 and Pr = on the left and 
right boundaries of Q, along with no flow (zero Neumann) conditions on the top 
and bottom boundaries. We also assume that there is no external forcing, i.e., that 
^ = 0. Table 1 shows the relative error of the dowscaled velocities obtained by 
the postprocessing procedure in Sect. 4. We observe that for the deterministic per- 
meability, an appropriate choice of L^ (the number of basis functions chosen for 
each respective node) significantly reduces the error associated with the down- 
scaled velocities. For the random field, we see that the error reduction is less 
pronounced, yet that a clear error decline is still evident. This diff'erence is due to 
the nature of the spectrum associated with the eigenvalue problem in Eq. (3.10) 
and, in particular, the eigenvalue behavior corresponding to the deterministic field 
lends itself to a more pronounced error decline (see, e.g., [10]). For this initial set 
of results all MsFEM/GMsFEM solutions were computed on 10 x 10 coarse mesh. 
For a visual representation of the improvement, the computed velocity is plot- 
ted on global and localized regions for the deterministic and random permeability 
fields in Figs. 5 and 6, respectively. In either figure, the left plots show the ref- 
erence fine scale velocity on the global domain, and on a subregion where the 
permeability has a large change in value. The center plots show the resulting 
downscaled velocity on the same regions resulting from the MsFEM (i.e., the 
case when L^ = 1). We note that significant diff'erences are noticeable between 
the reference values and this initial case. Finally, the right plots show the com- 
puted velocity corresponding GMsFEM with L^ = 4 on the same regions. Figs. 5 
and 6 clearly illustrate how the enrichment produces a more accurate represen- 
tation of the velocity on the fine scale. In addition, this increase in accuracy is 
particularly evident in regions where the permeability contrast is most extreme. 
We also note that there are some negative values for the horizontal velocity in the 
random permeability due to the channelized nature of the permeability in some 
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regions. 
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Figure 5: Velocity computed using three methods. The top two plots show the velocity profile on 
the whole domain using deterministic permeability with the reference on the left, L^ = 1 in the 
center, and L^ = 4 on the right. The bottom set of plots shows the velocity on the region shown in 
black and yellow on the top two plots. 



52. Two-phase convergence 

In solving Eq. (2.2) we use quadratic relative permeability curves kyy^ - S^ 
and kro = (1 - 5*)^, along withyu^; = 1 andyu^ = 5 for the fluid viscosities. For the 
initial condtion, the value at the left edge is set as S = 1 and we assume S (x, 0) = 
elsewhere. Results for application of the method to two-phase flow are shown in 
Figs. 7 and 8 using the same permeabilities from Fig. 4. Each of the plots shows 
the reference saturation in the first row at three time levels, L^ = 1 on the second 
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row at the same three time levels, L^ = 2 in the third row, and L^ = 4 in the last 
row. The improvement in the deterministic case is significant as can be seen in 
the Fig. 7 and verified with the relative error shown in Fig. 9. And as expected, 
for the random field we see a noticeable (but less pronounced) error decline from 
the profiles in Figs. 8, and the relative error plots in Fig. 9. We mention that in 
Figs. 7 and 8, the solutions corresponding to the case when L^ = 4 are essentially 
indistinguishable from the fine scale reference solutions. We reiterate that the size 
of the resulting linear systems from the pressure equation in Eq. (2.1) will depend 
on L^ and the coarse mesh size. The system for the reference cases are of size 
10201^ for the deterministic permeability and 14641^ for the random permeability. 
Since the coarse mesh is 10 x 10 for both cases, using L^ = 1, L^ = 2, and L^ = 4 
results in linear systems of size 121^, 202^, and 364^, respectively. 

The L^ error of the saturation profiles offered in Figs. 7 and 8 is presented 
in Fig. 9. As was evident from the velocity results, the error is improved by 
increasing the number of enrichment functions up to a certain threshhold, and 
then the reduction is minimal as more functions are added. This suggests the 
number of enrichment functions should be chosen to minimize the overhead in 
the calculation. As seen above, the size of the linear system to be solved grows 
quickly as the number of functions is increased and, at some point, there is little 
to be gain by adding to the enrichment. 

For a comparison on the mesh size, the deterministic permeability was used to 
simulate two-phase flow with L^ = 4 using a 10 x 10 and 25 x 25 coarse mesh. 
Analogous results are presented for the random field, except that a 30 x 30 coarse 
mesh was used for the refinement. The results are shown in Fig. 10. As expected, 
the error is reduced with a finer coarse mesh at the same level of enrichment. This 
is due to the fact that the coarse mesh is finer and the behavior of the permeability 
is captured to an acceptable degree by fewer basis functions. 

6. Concluding remarks 

In this paper we present a method for the construction of locally conservative 
flux fields from GMsFEM pressure solutions. The method hinges on an element- 
based postprocessing procedure in which an independent set of 4 x 4 linear sys- 
tems need to be solved in order to compute the flux values. In order to test the 
performance of the method we consider two permeability coefficients that exhibit 
distinct classes of heterogeneity. In addition, we apply the proposed method a 
two-phase flow model in which the flux values are coupled to a hyperbolic trans- 
port equation. The increase in accuracy associated with the computation of the 
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GMsFEM pressure solutions is inherited by the associated flux fields and satura- 
tion solutions, and is shown to closely depend on the size of the reduced-order 
systems. In particular, the addition of more basis functions to the enriched, mul- 
tiscale solution space yields solutions that more accurately capture the behavior 
of the fine scale model. In the future we wish to address related techniques in 
which the computation of conservative flux fields is built within the framework of 
generlized multiscale finite volume element methods. 

References 

References 

[1] J. Aarnes, S. Krogstad, and K. Lie, A hierarchical multiscale method for 
two-phase flow based on upon mixed finite elements and nonuniform coarse 
grids, SIAM Multiscale Model. Sim., 5 (2006), pp. 337-363. 

[2] T. Arbogast, G. Pencheva, M. Wheeler, and I. Yotov, A multiscale mor- 
tar mixed finite element method, SIAM Multiscale Model. Sim., 6 (2007), 
pp. 319-346. 

[3] K. Aziz and A. Settari, Petroleum reservoir simulation. Applied Science 
Pubhshers, 1979. 

[4] L. Bush and V. Ginting, On the application of the continuous galerkin finite 
element method for conservation problems (submitted), (2013). 

[5] J. Chu, Y. Efendiev, V. Ginting, and T. Hou, Flow based oversampling tech- 
nique for multiscale finite element methods. Advances in Water Resources, 
31 (2008), pp. 599 -608. 

[6] B. CocKBURN, J. GoPALAKRiSHNAN, AND H. Wang, Locally Conservative fluxcs 
for the continuous galerkin method, SIAM J. Numer. Anal., 45 (2007), 
pp. 1742-1776. 

[7] M. Dryja, On discontinuous galerkin methods for elliptic problems with dis- 
continuous coefficients, Comput. Methods Appl. Math., 3 (2003), pp. 76-85. 

[8] L. DuRLOFSKY, Y. Efendiev, and V. Ginting, An adaptive localglobal mul- 
tiscale finite volume element method for two-phase flow simulations. Ad- 
vances in Water Resources, 30 (2007), pp. 576 - 588. 



21 



[9] Y. Efendiev, J. Galvis, G. Li, and M. Presho, Generalized multiscale finite 
element methods, nonlinear elliptic equations (submitted), (2013). 

[10] Y. Efendiev, J. Galvis, and X. Wu, Multiscale finite element methods for 
high-contrast problems using local spectral basis functions, J. Comput. 
Phys., 230 (2011), pp. 937-955. 

[11] Y Efendiev, V. Ginting, T. Hou, and R. Ewing, Accurate multiscale finite ele- 
ment methods for two-phase flow simulations, J. Comput. Phys., 220 (2006), 
pp. 155-174. 

[12] Y Efendiev and T. Hou, Multiscale Finite Element Methods, Theory and 
Applications, Springer, New York, 2009. 

[13] L. C. Evans, Partial differential equations, vol. 19 of Graduate Studies in 
Mathematics, American Mathematical Society, Providence, RI, second ed., 
2010. 

[14] J. Galvis and Y Efendiev, Domain decomposition preconditioners for multi- 
scale flows in high contrast media, SIAM Multiscale Model. Sim., 8 (2010), 
pp. 1461-1483. 

[15] , Domain decomposition preconditioners for multiscale flows in high 

contrast media, reduced dimension coarse spaces, SIAM Multiscale Model. 
Sim., 8 (2010), pp. 1621-1644. 

[16] T. Hou AND X. Wu, A multiscale finite element method for elliptic problems 
in composite materials and porous media, J. Comput. Phys., 134 (1997), 
pp. 169-189. 

[17] T. Y Hou AND X.-H. Wu, A multiscale finite element method for elliptic 
problems in composite materials and porous media, J. Comput. Phys., 134 
(1997), pp. 169-189. 

[18] T. Y Hou, X.-H. Wu, and Z. Cai, Convergence of a multiscale finite ele- 
ment method for elliptic problems with rapidly oscillating coefficients. Math. 
Comp., 68 (1999), pp. 913-943. 

[19] T. Hughes, G. Engel, L. Mazzei, and M. Larson, The continuous galerkin 
method is locally conservative, J. Comput. Phys., 163 (2000), pp. 467-488. 



22 



[20] T. Hughes, G. Feijoo, L. Mazzei, and J. Quincy, The variational multiscale 
method - a paradigm for computational mechanics, Comput. Methods Appl. 
Mech. Engrg., 166 (1998), pp. 3-24. 

[21] P. Jenny, S. Lee, and H. Tchelepi, Multi-scale finite volume method for el- 
liptic problems in subsurface flow simulation, J. Comput. Phys., 187 (2003), 
pp. 47-67. 

[22] A. LouLA, F. RocHiNHA, AND M. MuRAD, Highcr-ordcr gradient post- 
processings for second-order elliptic problems, Comput. Methods Appl. 
Mech. Engrg., 128 (1995), pp. 361-381. 

[23] P.-A. Raviart and J. M. Thomas, A mixed finite element method for 2nd 
order elliptic problems, in Mathematical aspects of finite element meth- 
ods (Proc. Conf., Consigho Naz. delle Ricerche (C.N.R.), Rome, 1975), 
Springer, Berlin, 1977, pp. 292-315. Lecture Notes in Math., Vol. 606. 

[24] S. Sun and J. Liu, A locally conservative finite element method based on 
piecewise constant enrichment of the continuous galerkin method, SIAM J. 
Sci. Comput, 31 (2009), pp. 2528-2548. 

[25] J. W. Thomas, Numerical partial differential equations, vol. 33 of Texts in 
Applied Mathematics, Springer- Verlag, New York, 1999. Conservation laws 
and elliptic equations. 



23 




r''/^\ 



\ * 




'\ 



M\ 



y 



240 

140 

40 

-60 

-160 

-260 

-360 





^H 








K 


■H 








iH 


^^1 








HH 








^M 


^^1 


1 |ga 


§ 


t ^ 


sRh 


■■ 










R 


9 




1 


■Sn 


H 


u 








,a< 


■■■Uiii 


1 


^H 




u 


idhi 




■ 


H 


iH 




■ 


^^^^Hl 


n 


■ 




11 


VHSa 










■■■ 


J|B9| 


:,,- 


BB 




IT 





Figure 6: Velocity computed using three methods. The top two plots show the velocity profile on 
the whole domain using random permeability with the reference on the left, L^ = 1 in the center, 
and L^ = 4 on the right. The bottom set of plots shows the velocity on the region shown in black 
and yellow on the top two plots. 
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Figure 7: Deterministic permeability. The reference saturation is shown in the top row at three 
different time levels. The second through fourth row are L^ = 1, L^ = 2, and L^ = 4 respectively. 
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Figure 8: Random permeability. The reference saturation is shown in the top row at three different 
time levels. The second through fourth row are L^ = 1, L^ = 2, and L^ = 4 respectively. 
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Figure 9: Comparison of the L^ error of the saturation for deterministic (left) and random (right) 
as a function of time. 
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Figure 10: Error comparison as the coarse mesh is changed with L^ = 4 and both permeabilities. 
Filled circles show the error behavior as a function of time with a 10 x 10 coarse mesh. Squares 
show the error behavior with a finer coarse mesh. 
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